knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 2000 # 10'000 per chain to achieve 40'000
warmup = 1000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = as.character(iterations)
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]

Constructing scales reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,24)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,24+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,47)
  )

rows_to_pack_hu_ordinal <- list(
  "Intercepts" = c(1,7),
  "Conditional Within-Person Effects" = c(3+5,10+5),
  "Conditional Between-Person Effects" = c(11+5,17+5),
  
  "Hurdle Within-Person Effects" = c(18+5,25+5),
  "Hurdle Between-Person Effects" = c(26+5,32+5),
  
  "Random Effects" = c(33+5, 46+5), 
  "Additional Parameters" = c(47+5,47+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

WITHOUT PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_Nopushing", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -10536.2 166.2
## p_loo       144.3   4.8
## looic     21072.3 332.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  512     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 8, observations = 3736, p-value = 0.46
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000267666 0.003211991
## sample estimates:
## outlier frequency (expected: 0.00152837259100642 ) 
##                                        0.002141328
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 48.72*** 3.08 [42.92, 55.10] 1.000 [0.92, 1.08] 0.000 1.002 964 1749
Hurdle Intercept 0.82 0.14 [ 0.59, 1.13] 0.882 [0.84, 1.20] 0.433 1.005 552 1254
Conditional Within-Person Effects
Daily persuasion experienced 1.03 0.03 [ 0.97, 1.09] 0.853 [0.92, 1.08] 0.965 1.001 1471 2332
Daily persuasion utilized (partner’s view) 1.03 0.02 [ 0.98, 1.08] 0.877 [0.92, 1.08] 0.987 1.002 2066 2263
Daily pressure experienced 0.90* 0.04 [ 0.81, 0.99] 0.983 [0.92, 1.08] 0.260 1.000 3707 2336
Daily pressure utilized (partner’s view) 0.94 0.04 [ 0.86, 1.02] 0.928 [0.92, 1.08] 0.634 1.002 4371 3030
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 0.98 0.06 [ 0.87, 1.11] 0.623 [0.92, 1.08] 0.766 1.000 5745 2727
Daily weartime NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.11 0.15 [ 0.85, 1.47] 0.775 [0.92, 1.08] 0.336 1.004 817 1463
Mean persuasion utilized (partner’s view) 1.09 0.15 [ 0.83, 1.45] 0.737 [0.92, 1.08] 0.354 1.007 812 1214
Mean pressure experienced 1.12 0.18 [ 0.81, 1.55] 0.752 [0.92, 1.08] 0.287 1.005 1015 1888
Mean pressure utilized (partner’s view) 0.91 0.16 [ 0.64, 1.28] 0.708 [0.92, 1.08] 0.305 1.003 1145 1670
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.61*** 0.10 [ 1.43, 1.84] 1.000 [0.84, 1.20] 0.000 1.001 1969 2771
Hu Daily persuasion utilized (partner’s view) 1.43*** 0.08 [ 1.28, 1.61] 1.000 [0.84, 1.20] 0.000 1.001 3179 3039
Hu Daily pressure experienced 1.01 0.15 [ 0.76, 1.39] 0.525 [0.84, 1.20] 0.779 1.000 3515 2678
Hu Daily pressure utilized (partner’s view) 1.74*** 0.32 [ 1.26, 2.83] 1.000 [0.84, 1.20] 0.009 1.001 2286 1405
Hu Daily pushing experienced NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Day 0.92 0.12 [ 0.71, 1.20] 0.730 [0.84, 1.20] 0.759 1.001 7057 2959
Hu Daily weartime NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.56 0.55 [ 0.76, 3.18] 0.889 [0.84, 1.20] 0.183 1.011 579 1175
Hu Mean persuasion utilized (partner’s view) 1.55 0.55 [ 0.76, 3.19] 0.881 [0.84, 1.20] 0.188 1.009 606 1142
Hu Mean pressure experienced 0.33** 0.14 [ 0.15, 0.77] 0.996 [0.84, 1.20] 0.013 1.007 782 1578
Hu Mean pressure utilized (partner’s view) 0.58 0.24 [ 0.25, 1.33] 0.911 [0.84, 1.20] 0.139 1.007 803 1296
Hu Mean pushing experienced NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.32 0.04 [0.24, 0.42] NA NA NA 1.008 985 1765
sd(Hurdle Intercept) 0.89 0.12 [0.70, 1.18] NA NA NA 1.001 902 1729
sd(Daily persuasion experienced) 0.12 0.02 [0.08, 0.17] NA NA NA 1.003 1813 2949
sd(Daily persuasion utilized (partner’s view)) 0.09 0.02 [0.05, 0.14] NA NA NA 1.006 2186 2589
sd(Daily pressure experienced) 0.08 0.07 [0.00, 0.25] NA NA NA 1.001 1377 1876
sd(Daily pressure utilized (partner’s view)) 0.07 0.05 [0.00, 0.19] NA NA NA 1.003 1648 1853
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion experienced) 0.21 0.07 [0.08, 0.37] NA NA NA 1.002 1341 977
sd(Hu Daily persuasion utilized (partner’s view)) 0.19 0.07 [0.05, 0.34] NA NA NA 1.001 1429 1309
sd(Hu Daily pressure experienced) 0.23 0.20 [0.01, 0.80] NA NA NA 1.003 1209 1877
sd(Hu Daily pressure utilized (partner’s view)) 0.31 0.25 [0.02, 1.01] NA NA NA 1.004 1222 1709
sd(Hu Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.69 0.01 [0.67, 0.71] NA NA NA 1.001 5060 3213
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.77), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.77), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.82), b_hu_pressure_self_cb and
##   b_hu_persuasion_self_cb (r = 0.81), b_hu_pressure_partner_cb and
##   b_hu_persuasion_partner_cb (r = 0.81). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 157 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 173 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.011

$persuasion_partner_cw

## Warning: Removed 143 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 145 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00867

$pressure_self_cw

## Warning: Removed 46 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 55 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.016

$pressure_partner_cw

## Warning: Removed 49 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 85 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0275

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_nopushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -18783.3  68.4
## p_loo        74.9   3.4
## looic     37566.6 136.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.7]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 27, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001498352 0.005394067
## sample estimates:
## outlier frequency (expected: 0.00329937069223854 ) 
##                                          0.0080911
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 116.98*** 6.11 [105.53, 130.15] 1.000 [0.94, 1.07] 0.000 1.003 583 1345
Within-Person Effects
Daily persuasion experienced 1.04* 0.01 [ 1.01, 1.07] 0.989 [0.94, 1.07] 0.975 1.001 2606 2982
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.99, 1.06] 0.936 [0.94, 1.07] 0.995 1.000 3058 3179
Daily pressure experienced 0.95 0.03 [ 0.89, 1.02] 0.935 [0.94, 1.07] 0.678 1.000 4455 3085
Daily pressure utilized (partner’s view) 0.99 0.03 [ 0.93, 1.06] 0.627 [0.94, 1.07] 0.936 1.001 4236 2961
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 0.97 0.03 [ 0.91, 1.04] 0.774 [0.94, 1.07] 0.856 1.001 7048 2972
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 1.001 3742 2490
Between-Person Effects
Mean persuasion experienced 1.11 0.14 [ 0.85, 1.42] 0.788 [0.94, 1.07] 0.279 1.007 529 1039
Mean persuasion utilized (partner’s view) 1.04 0.13 [ 0.79, 1.34] 0.611 [0.94, 1.07] 0.374 1.008 538 1075
Mean pressure experienced 0.89 0.13 [ 0.68, 1.20] 0.782 [0.94, 1.07] 0.244 1.004 669 1436
Mean pressure utilized (partner’s view) 1.02 0.14 [ 0.78, 1.37] 0.578 [0.94, 1.07] 0.362 1.005 643 1395
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.912 [0.94, 1.07] 1.000 1.000 4781 3605
Random Effects
sd(Intercept) 0.30 0.04 [0.23, 0.39] NA NA NA 1.005 957 1591
sd(Daily persuasion experienced) 0.05 0.01 [0.03, 0.08] NA NA NA 1.002 2335 2032
sd(Daily persuasion utilized (partner’s view)) 0.06 0.02 [0.03, 0.09] NA NA NA 1.001 1909 2207
sd(Daily pressure experienced) 0.05 0.04 [0.00, 0.14] NA NA NA 1.001 1627 2057
sd(Daily pressure utilized (partner’s view)) 0.04 0.03 [0.00, 0.12] NA NA NA 1.000 2070 1630
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.58 0.01 [0.56, 0.59] NA NA NA 1.004 7025 2589
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.88), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.85), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.77), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.87). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_nopushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
check_brms(mood_gauss, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5189.1  59.2
## p_loo        68.2   3.0
## looic     10378.3 118.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  750     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 32, observations = 3736, p-value =
## 2.24e-11
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005865832 0.012070325
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.00856531
summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 3.70*** 0.11 [ 3.49, 3.90] 1.000 [-0.11, 0.11] 0.000 1.006 417 729
Within-Person Effects
Daily persuasion experienced 0.01 0.02 [-0.03, 0.05] 0.663 [-0.11, 0.11] 1.000 1.001 3672 2684
Daily persuasion utilized (partner’s view) 0.04 0.02 [-0.01, 0.08] 0.938 [-0.11, 0.11] 0.999 1.000 3341 3254
Daily pressure experienced -0.03 0.05 [-0.14, 0.07] 0.743 [-0.11, 0.11] 0.940 1.000 4889 2934
Daily pressure utilized (partner’s view) 0.01 0.05 [-0.10, 0.11] 0.559 [-0.11, 0.11] 0.965 1.001 3761 2704
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 0.25*** 0.06 [ 0.15, 0.36] 1.000 [-0.11, 0.11] 0.006 1.001 7106 2660
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.42 0.25 [-0.10, 0.91] 0.942 [-0.11, 0.11] 0.099 1.010 499 775
Mean persuasion utilized (partner’s view) 0.33 0.25 [-0.18, 0.82] 0.902 [-0.11, 0.11] 0.152 1.009 487 558
Mean pressure experienced -0.36 0.26 [-0.91, 0.19] 0.903 [-0.11, 0.11] 0.139 1.007 571 954
Mean pressure utilized (partner’s view) -0.29 0.27 [-0.82, 0.27] 0.851 [-0.11, 0.11] 0.180 1.007 549 835
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.07 [0.48, 0.79] NA NA NA 1.009 812 1382
sd(Daily persuasion experienced) 0.05 0.03 [0.00, 0.11] NA NA NA 1.005 1000 1287
sd(Daily persuasion utilized (partner’s view)) 0.08 0.03 [0.03, 0.14] NA NA NA 1.003 1363 1122
sd(Daily pressure experienced) 0.06 0.06 [0.00, 0.23] NA NA NA 1.001 1944 2092
sd(Daily pressure utilized (partner’s view)) 0.07 0.06 [0.00, 0.25] NA NA NA 1.002 2091 2255
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.96 0.01 [0.94, 0.98] NA NA NA 1.000 8235 2720
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.86), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.86), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_nopushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
check_brms(reactance_ordinal)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -683.2 31.6
## p_loo        59.6  4.5
## looic      1366.3 63.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     751   99.3%   224     
##    (0.7, 1]   (bad)        5    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 0.36
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002645503
## sample estimates:
## outlier frequency (expected: 0.000304232804232804 ) 
##                                         0.001322751
summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA
Intercept[1] 3.48*** 0.86 [ 2.16, 5.75] 1.000 [0.84, 1.20] 0.000 1.000 3268 2969
Intercept[2] 7.32*** 1.90 [ 4.48, 12.47] 1.000 [0.84, 1.20] 0.000 1.001 3373 3061
Intercept[3] 19.78*** 5.57 [ 11.49, 35.14] 1.000 [0.84, 1.20] 0.000 1.001 3544 3388
Intercept[4] 84.24*** 29.47 [ 45.16, 173.12] 1.000 [0.84, 1.20] 0.000 1.000 4264 3330
Intercept[5] 2760.30*** 1715.32 [929.55, 10049.83] 1.000 [0.84, 1.20] 0.000 1.000 5669 2968
Within-Person Effects
Daily persuasion experienced 0.85* 0.07 [ 0.72, 0.99] 0.979 [0.84, 1.20] 0.596 1.000 4158 3111
Daily persuasion utilized (partner’s view) 1.01 0.09 [ 0.84, 1.20] 0.548 [0.84, 1.20] 0.950 1.001 4250 2828
Daily pressure experienced 1.98** 0.34 [ 1.31, 2.76] 0.999 [0.84, 1.20] 0.011 1.001 2901 3115
Daily pressure utilized (partner’s view) 1.15 0.23 [ 0.72, 1.85] 0.749 [0.84, 1.20] 0.509 1.001 3088 2320
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 1.63 0.55 [ 0.87, 3.10] 0.929 [0.84, 1.20] 0.161 1.000 6064 2790
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.08 0.50 [ 0.42, 2.76] 0.567 [0.84, 1.20] 0.309 1.005 1276 2066
Mean persuasion utilized (partner’s view) 0.79 0.42 [ 0.29, 2.26] 0.665 [0.84, 1.20] 0.243 1.004 1346 2146
Mean pressure experienced 3.30* 1.70 [ 1.21, 9.23] 0.990 [0.84, 1.20] 0.019 1.004 1610 2466
Mean pressure utilized (partner’s view) 1.23 0.70 [ 0.41, 3.54] 0.637 [0.84, 1.20] 0.226 1.003 1666 2712
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.78 0.19 [0.46, 1.22] NA NA NA 1.001 1742 2619
sd(Daily persuasion experienced) 0.16 0.12 [0.01, 0.44] NA NA NA 1.008 749 1628
sd(Daily persuasion utilized (partner’s view)) 0.18 0.13 [0.01, 0.45] NA NA NA 1.002 1314 1689
sd(Daily pressure experienced) 0.47 0.22 [0.07, 1.02] NA NA NA 1.007 1083 919
sd(Daily pressure utilized (partner’s view)) 0.31 0.31 [0.02, 1.37] NA NA NA 1.002 1028 1907
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.83),
##   b_pressure_self_cb and b_persuasion_self_cb (r = 0.79),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.84). This might
##   lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="03_SensitivityPushingSeparate_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_nopushing_", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
check_brms(is_reactance)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 31 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -361.7 15.6
## p_loo        63.7  4.9
## looic       723.3 31.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     725   95.9%   78      
##    (0.7, 1]   (bad)       30    4.0%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.6663
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0003205434 0.0095234999
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002645503
summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 0.32*** 0.08 [0.18, 0.53] 1.000 [0.83, 1.20] 0.000 1.002 2308 2705
Within-Person Effects
Daily persuasion experienced 0.86 0.08 [0.71, 1.01] 0.964 [0.83, 1.20] 0.612 1.000 2887 2614
Daily persuasion utilized (partner’s view) 1.10 0.14 [0.86, 1.45] 0.777 [0.83, 1.20] 0.729 1.001 2245 2471
Daily pressure experienced 2.09* 0.58 [1.19, 4.20] 0.991 [0.83, 1.20] 0.024 1.001 1589 1427
Daily pressure utilized (partner’s view) 1.34 0.53 [0.57, 3.87] 0.780 [0.83, 1.20] 0.276 1.001 1658 1539
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 1.86 0.67 [0.89, 3.88] 0.955 [0.83, 1.20] 0.092 1.001 4640 2592
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.86 1.07 [0.58, 5.93] 0.859 [0.83, 1.20] 0.144 1.005 1143 2107
Mean persuasion utilized (partner’s view) 1.03 0.66 [0.31, 3.66] 0.520 [0.83, 1.20] 0.233 1.004 1265 2190
Mean pressure experienced 17.76*** 16.40 [3.37, 123.01] 1.000 [0.83, 1.20] 0.000 1.002 1423 2478
Mean pressure utilized (partner’s view) 1.21 1.21 [0.17, 7.77] 0.574 [0.83, 1.20] 0.136 1.001 1712 2485
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.21 0.25 [0.80, 1.78] NA NA NA 1.001 1396 1835
sd(Daily persuasion experienced) 0.19 0.13 [0.01, 0.48] NA NA NA 1.004 958 1649
sd(Daily persuasion utilized (partner’s view)) 0.43 0.17 [0.12, 0.83] NA NA NA 1.003 1279 1064
sd(Daily pressure experienced) 0.79 0.53 [0.07, 2.12] NA NA NA 1.004 481 715
sd(Daily pressure utilized (partner’s view)) 0.74 0.64 [0.04, 2.71] NA NA NA 1.002 1030 1680
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsNoPushing.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 48.72*** [42.92, 55.10] 1.000 116.98*** [105.53, 130.15] 1.000 3.70*** [ 3.49, 3.90] 1.000 NA NA NA 0.32*** [0.18, 0.53] 1.000
Hurdle Intercept 0.82 [ 0.59, 1.13] 0.882 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.03 [ 0.97, 1.09] 0.853 1.04* [ 1.01, 1.07] 0.989 0.01 [-0.03, 0.05] 0.663 0.85* [ 0.72, 0.99] 0.979 0.86 [0.71, 1.01] 0.964
Daily persuasion utilized (partner’s view) 1.03 [ 0.98, 1.08] 0.877 1.02 [ 0.99, 1.06] 0.936 0.04 [-0.01, 0.08] 0.938 1.01 [ 0.84, 1.20] 0.548 1.10 [0.86, 1.45] 0.777
Daily pressure experienced 0.90* [ 0.81, 0.99] 0.983 0.95 [ 0.89, 1.02] 0.935 -0.03 [-0.14, 0.07] 0.743 1.98** [ 1.31, 2.76] 0.999 2.09* [1.19, 4.20] 0.991
Daily pressure utilized (partner’s view) 0.94 [ 0.86, 1.02] 0.928 0.99 [ 0.93, 1.06] 0.627 0.01 [-0.10, 0.11] 0.559 1.15 [ 0.72, 1.85] 0.749 1.34 [0.57, 3.87] 0.780
Daily pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Day 0.98 [ 0.87, 1.11] 0.623 0.97 [ 0.91, 1.04] 0.774 0.25*** [ 0.15, 0.36] 1.000 1.63 [ 0.87, 3.10] 0.929 1.86 [0.89, 3.88] 0.955
Daily weartime NA NA NA 1.00*** [ 1.00, 1.00] 1.000 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.11 [ 0.85, 1.47] 0.775 1.11 [ 0.85, 1.42] 0.788 0.42 [-0.10, 0.91] 0.942 1.08 [ 0.42, 2.76] 0.567 1.86 [0.58, 5.93] 0.859
Mean persuasion utilized (partner’s view) 1.09 [ 0.83, 1.45] 0.737 1.04 [ 0.79, 1.34] 0.611 0.33 [-0.18, 0.82] 0.902 0.79 [ 0.29, 2.26] 0.665 1.03 [0.31, 3.66] 0.520
Mean pressure experienced 1.12 [ 0.81, 1.55] 0.752 0.89 [ 0.68, 1.20] 0.782 -0.36 [-0.91, 0.19] 0.903 3.30* [ 1.21, 9.23] 0.990 17.76*** [3.37, 123.01] 1.000
Mean pressure utilized (partner’s view) 0.91 [ 0.64, 1.28] 0.708 1.02 [ 0.78, 1.37] 0.578 -0.29 [-0.82, 0.27] 0.851 1.23 [ 0.41, 3.54] 0.637 1.21 [0.17, 7.77] 0.574
Mean pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.912 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.61*** [ 1.43, 1.84] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 1.43*** [ 1.28, 1.61] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 1.01 [ 0.76, 1.39] 0.525 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 1.74*** [ 1.26, 2.83] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.92 [ 0.71, 1.20] 0.730 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.56 [ 0.76, 3.18] 0.889 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 1.55 [ 0.76, 3.19] 0.881 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 0.33** [ 0.15, 0.77] 0.996 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.58 [ 0.25, 1.33] 0.911 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.32 [0.24, 0.42] NA 0.30 [0.23, 0.39] NA 0.61 [0.48, 0.79] NA 0.78 [0.46, 1.22] NA 1.21 [0.80, 1.78] NA
sd(Hurdle Intercept) 0.89 [0.70, 1.18] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.12 [0.08, 0.17] NA 0.05 [0.03, 0.08] NA 0.05 [0.00, 0.11] NA 0.16 [0.01, 0.44] NA 0.19 [0.01, 0.48] NA
sd(Daily persuasion utilized (partner’s view)) 0.09 [0.05, 0.14] NA 0.06 [0.03, 0.09] NA 0.08 [0.03, 0.14] NA 0.18 [0.01, 0.45] NA 0.43 [0.12, 0.83] NA
sd(Daily pressure experienced) 0.08 [0.00, 0.25] NA 0.05 [0.00, 0.14] NA 0.06 [0.00, 0.23] NA 0.47 [0.07, 1.02] NA 0.79 [0.07, 2.12] NA
sd(Daily pressure utilized (partner’s view)) 0.07 [0.00, 0.19] NA 0.04 [0.00, 0.12] NA 0.07 [0.00, 0.25] NA 0.31 [0.02, 1.37] NA 0.74 [0.04, 2.71] NA
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion experienced) 0.21 [0.08, 0.37] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.19 [0.05, 0.34] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.23 [0.01, 0.80] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.31 [0.02, 1.01] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.69 [0.67, 0.71] NA 0.58 [0.56, 0.59] NA 0.96 [0.94, 0.98] NA NA NA NA NA NA NA

ONLY PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_onlypushing", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
## Start sampling
pa_sub_digest <- digest::digest(pa_sub)
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 1342 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5692.9  78.2
## p_loo        97.8   4.1
## looic     11385.8 156.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1341  99.9%   633     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 4, observations = 1342, p-value = 0.26
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003371833
## sample estimates:
## outlier frequency (expected: 0.00149776453055142 ) 
##                                        0.002980626
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 57.28*** 4.53 [49.10, 66.97] 1.000 [0.93, 1.08] 0.000 1.002 1327 2040
Hurdle Intercept 3.79*** 0.98 [ 2.27, 6.26] 1.000 [0.84, 1.20] 0.000 1.005 1088 2463
Conditional Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 0.97 0.03 [ 0.91, 1.03] 0.856 [0.93, 1.08] 0.931 1.002 4427 2662
Daily pushing utilized (partner’s view) 0.97 0.03 [ 0.90, 1.03] 0.847 [0.93, 1.08] 0.907 1.001 3679 2873
Day 0.91 0.07 [ 0.79, 1.06] 0.894 [0.93, 1.08] 0.398 1.000 8744 2600
Daily weartime NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 1.32* 0.18 [ 1.00, 1.71] 0.977 [0.93, 1.08] 0.060 1.005 1217 2176
Mean pushing utilized (partner’s view) 1.18 0.15 [ 0.90, 1.53] 0.883 [0.93, 1.08] 0.222 1.005 1265 2076
Mean weartime NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 1.55** 0.23 [ 1.17, 2.24] 0.998 [0.84, 1.20] 0.035 1.002 3127 2531
Hu Daily pushing utilized (partner’s view) 1.84*** 0.30 [ 1.38, 2.68] 1.000 [0.84, 1.20] 0.002 1.001 2766 2989
Hu Day 0.61 0.15 [ 0.37, 1.00] 0.974 [0.84, 1.20] 0.100 1.003 7632 2768
Hu Daily weartime NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 0.46** 0.14 [ 0.26, 0.81] 0.997 [0.84, 1.20] 0.018 1.001 2675 2824
Hu Mean pushing utilized (partner’s view) 0.62 0.18 [ 0.33, 1.08] 0.952 [0.84, 1.20] 0.136 1.001 2568 2699
Hu Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.37 0.06 [0.29, 0.50] NA NA NA 1.002 1283 2370
sd(Hurdle Intercept) 1.25 0.19 [0.93, 1.74] NA NA NA 1.002 1448 2426
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.06 0.04 [0.00, 0.16] NA NA NA 1.001 870 1486
sd(Daily pushing utilized (partner’s view)) 0.09 0.04 [0.01, 0.18] NA NA NA 1.004 977 1142
sd(Hu Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.46 0.24 [0.07, 1.02] NA NA NA 1.000 947 1689
sd(Hu Daily pushing utilized (partner’s view)) 0.51 0.24 [0.10, 1.13] NA NA NA 1.001 1088 1582
Additional Parameters
sigma 0.67 0.02 [0.64, 0.71] NA NA NA 1.002 6374 2773
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$pushing_self_cw

## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0172

$pushing_partner_cw

## Warning: Removed 130 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 180 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.022

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
   
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_onlypushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 1214 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -6914.2 39.4
## p_loo        58.6  4.3
## looic     13828.4 78.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 8, observations = 1214, p-value = 0.16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0008237232 0.0074135091
## sample estimates:
## outlier frequency (expected: 0.00339373970345964 ) 
##                                        0.006589786
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 124.24*** 7.52 [109.69, 139.58] 1.000 [0.94, 1.06] 0.000 1.002 967 1634
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 1.03 0.03 [ 0.97, 1.08] 0.807 [0.94, 1.06] 0.914 1.001 3221 2733
Daily pushing utilized (partner’s view) 1.02 0.02 [ 0.97, 1.06] 0.797 [0.94, 1.06] 0.979 1.000 5065 2748
Day 0.97 0.05 [ 0.87, 1.08] 0.718 [0.94, 1.06] 0.654 1.000 7350 3336
Daily weartime 1.00** 0.00 [ 1.00, 1.00] 0.998 [0.94, 1.06] 1.000 1.001 4680 2898
Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 0.97 0.06 [ 0.86, 1.10] 0.680 [0.94, 1.06] 0.635 1.001 2223 2645
Mean pushing utilized (partner’s view) 1.04 0.07 [ 0.92, 1.18] 0.739 [0.94, 1.06] 0.578 1.002 2058 2692
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.812 [0.94, 1.06] 1.000 1.000 3879 3246
Random Effects
sd(Intercept) 0.31 0.04 [0.24, 0.42] NA NA NA 1.005 890 1790
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.09 0.04 [0.02, 0.17] NA NA NA 1.004 887 821
sd(Daily pushing utilized (partner’s view)) 0.03 0.03 [0.00, 0.10] NA NA NA 1.001 1398 1731
Additional Parameters
sigma 0.54 0.01 [0.52, 0.57] NA NA NA 1.001 6040 3037
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_onlypushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
check_brms(mood_gauss, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 1342 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1798.8 39.9
## p_loo        50.3  3.0
## looic      3597.5 79.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 13, observations = 1342, p-value =
## 4.837e-06
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005167734 0.016508158
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.009687034
summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 3.78*** 0.12 [ 3.56, 4.00] 1.000 [-0.11, 0.11] 0.000 1.013 560 1111
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 0.02 0.04 [-0.05, 0.08] 0.670 [-0.11, 0.11] 0.995 1.002 2823 2591
Daily pushing utilized (partner’s view) 0.09* 0.04 [ 0.01, 0.17] 0.984 [-0.11, 0.11] 0.715 1.002 1597 2207
Day 0.28*** 0.09 [ 0.11, 0.46] 1.000 [-0.11, 0.11] 0.023 1.001 3871 3088
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 0.12 0.11 [-0.08, 0.32] 0.880 [-0.11, 0.11] 0.433 1.003 1269 1954
Mean pushing utilized (partner’s view) 0.08 0.11 [-0.12, 0.30] 0.786 [-0.11, 0.11] 0.554 1.002 1216 1958
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.08 [0.48, 0.79] NA NA NA 1.011 648 1448
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.07 0.04 [0.00, 0.16] NA NA NA 1.000 1409 1126
sd(Daily pushing utilized (partner’s view)) 0.14 0.04 [0.06, 0.24] NA NA NA 1.002 1329 828
Additional Parameters
sigma 0.91 0.02 [0.87, 0.94] NA NA NA 1.001 3887 2536
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_onlypushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
check_brms(reactance_ordinal)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -363.6 23.5
## p_loo        33.6  2.9
## looic       727.1 47.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.6]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 403, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.00248139
## sample estimates:
## outlier frequency (expected: 0.000297766749379653 ) 
##                                                   0
summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA
Intercept[1] 5.44*** 1.57 [ 3.06, 10.06] 1.000 [0.84, 1.20] 0.000 1.001 3521 2882
Intercept[2] 9.89*** 3.14 [ 5.45, 18.74] 1.000 [0.84, 1.20] 0.000 1.002 3598 2951
Intercept[3] 23.10*** 7.89 [ 12.38, 46.79] 1.000 [0.84, 1.20] 0.000 1.001 3933 2899
Intercept[4] 81.88*** 35.09 [ 38.53, 199.63] 1.000 [0.84, 1.20] 0.000 1.000 4569 3054
Intercept[5] 771.25*** 625.85 [202.52, 4855.18] 1.000 [0.84, 1.20] 0.000 1.001 5246 2545
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 1.31* 0.14 [ 1.04, 1.62] 0.987 [0.84, 1.20] 0.214 1.001 3641 2884
Daily pushing utilized (partner’s view) 0.96 0.13 [ 0.75, 1.25] 0.615 [0.84, 1.20] 0.810 1.002 4044 2794
Day 1.60 0.73 [ 0.68, 3.95] 0.854 [0.84, 1.20] 0.182 1.002 4822 3245
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 2.32** 0.65 [ 1.35, 3.99] 1.000 [0.84, 1.20] 0.005 1.001 3394 2987
Mean pushing utilized (partner’s view) 1.25 0.35 [ 0.70, 2.16] 0.783 [0.84, 1.20] 0.365 1.001 3738 2562
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.85 0.21 [0.49, 1.34] NA NA NA 1.001 1569 2204
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.17 0.14 [0.01, 0.50] NA NA NA 1.004 1241 1651
sd(Daily pushing utilized (partner’s view)) 0.16 0.14 [0.01, 0.57] NA NA NA 1.003 1626 1714
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.82). This might lead
##   to inappropriate results. See 'Details' in '?rope'.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_onlypushing", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
check_brms(is_reactance)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -216.1 11.0
## p_loo        39.3  3.0
## looic       432.3 22.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     400   99.3%   458     
##    (0.7, 1]   (bad)        3    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 403, p-value = 0.5534
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  6.282137e-05 1.374725e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.00248139
summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 0.40** 0.11 [0.23, 0.68] 1.000 [0.83, 1.20] 0.002 1.000 3095 2830
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 1.37* 0.18 [1.06, 1.87] 0.993 [0.83, 1.20] 0.154 1.000 3238 2410
Daily pushing utilized (partner’s view) 1.16 0.26 [0.80, 2.12] 0.769 [0.83, 1.20] 0.509 1.000 2310 2170
Day 1.78 0.89 [0.69, 4.81] 0.874 [0.83, 1.20] 0.160 1.000 6984 2927
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 5.39*** 2.58 [2.38, 14.94] 1.000 [0.83, 1.20] 0.000 1.001 2293 2628
Mean pushing utilized (partner’s view) 2.37 1.15 [0.97, 6.82] 0.972 [0.83, 1.20] 0.057 1.001 2556 2640
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.52 0.31 [1.03, 2.24] NA NA NA 1.001 1431 2494
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.18 0.15 [0.01, 0.56] NA NA NA 1.003 1256 1876
sd(Daily pushing utilized (partner’s view)) 0.42 0.27 [0.03, 1.15] NA NA NA 1.002 1292 1879
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_self_cw

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsOnlyPushing.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 57.28*** [49.10, 66.97] 1.000 124.24*** [109.69, 139.58] 1.000 3.78*** [ 3.56, 4.00] 1.000 NA NA NA 0.40** [0.23, 0.68] 1.000
Hurdle Intercept 3.79*** [ 2.27, 6.26] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pushing experienced 0.97 [ 0.91, 1.03] 0.856 1.03 [ 0.97, 1.08] 0.807 0.02 [-0.05, 0.08] 0.670 1.31* [ 1.04, 1.62] 0.987 1.37* [1.06, 1.87] 0.993
Daily pushing utilized (partner’s view) 0.97 [ 0.90, 1.03] 0.847 1.02 [ 0.97, 1.06] 0.797 0.09* [ 0.01, 0.17] 0.984 0.96 [ 0.75, 1.25] 0.615 1.16 [0.80, 2.12] 0.769
Day 0.91 [ 0.79, 1.06] 0.894 0.97 [ 0.87, 1.08] 0.718 0.28*** [ 0.11, 0.46] 1.000 1.60 [ 0.68, 3.95] 0.854 1.78 [0.69, 4.81] 0.874
Daily weartime NA NA NA 1.00** [ 1.00, 1.00] 0.998 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pushing experienced 1.32* [ 1.00, 1.71] 0.977 0.97 [ 0.86, 1.10] 0.680 0.12 [-0.08, 0.32] 0.880 2.32** [ 1.35, 3.99] 1.000 5.39*** [2.38, 14.94] 1.000
Mean pushing utilized (partner’s view) 1.18 [ 0.90, 1.53] 0.883 1.04 [ 0.92, 1.18] 0.739 0.08 [-0.12, 0.30] 0.786 1.25 [ 0.70, 2.16] 0.783 2.37 [0.97, 6.82] 0.972
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.812 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 1.55** [ 1.17, 2.24] 0.998 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 1.84*** [ 1.38, 2.68] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.61 [ 0.37, 1.00] 0.974 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 0.46** [ 0.26, 0.81] 0.997 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 0.62 [ 0.33, 1.08] 0.952 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.37 [0.29, 0.50] NA 0.31 [0.24, 0.42] NA 0.61 [0.48, 0.79] NA 0.85 [0.49, 1.34] NA 1.52 [1.03, 2.24] NA
sd(Hurdle Intercept) 1.25 [0.93, 1.74] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.06 [0.00, 0.16] NA 0.09 [0.02, 0.17] NA 0.07 [0.00, 0.16] NA 0.17 [0.01, 0.50] NA 0.18 [0.01, 0.56] NA
sd(Daily pushing utilized (partner’s view)) 0.09 [0.01, 0.18] NA 0.03 [0.00, 0.10] NA 0.14 [0.06, 0.24] NA 0.16 [0.01, 0.57] NA 0.42 [0.03, 1.15] NA
sd(Hu Daily persuasion experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.46 [0.07, 1.02] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.51 [0.10, 1.13] NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.67 [0.64, 0.71] NA 0.54 [0.52, 0.57] NA 0.91 [0.87, 0.94] NA NA NA NA NA NA NA
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • fs (version 1.6.5; Hester J et al., 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()